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Electronic levels and electrical response of periodic molecular structures from 
plane-wave orbital-dependent calculations 

Yanli Li^ and Ismaila Dabo^'Lj 

^Department of Physics, Xiamen University, Xiamen 361005, Republic of China 
^Umversite Pans-Est, CERMICS, Projet Micmac ENPC-INRIA, 
6-8 avenue Blaise Pascal, 77455 Marne-la-Vallee Cedex 2, France 

Plane-wave electronic-structure predictions based upon orbital-dependent density-functional the- 
ory (OD-DFT) approximations, such as hybrid density-functional methods and self-interaction 
density-functional corrections, are severely affected by computational inaccuracies in evaluating elec- 
tron interactions in the plane-wave representation. These errors arise from divergence singularities 
in the plane-wave summation of electrostatic and exchange interaction contributions. Auxiliary- 
function corrections are reciprocal-space countercharge corrections that cancel plane-wave singu- 
larities through the addition of an auxiliary function to the point-charge electrostatic kernel that 
enters into the expression of interaction terms. At variance with real-space countercharge corrections 
that are employed in the context of density- functional theory (DFT), reciprocal-space corrections 
are computationally inexpensive, making them suited to more demanding OD-DFT calculations. 
Nevertheless, there exists much freedom in the choice of auxiliary functions and various definitions 
result in different levels of performance in eliminating plane- wave inaccuracies. In this work, we 
derive exact point-charge auxiliary functions for the description of molecular structures of arbitrary 
translational symmetry, including the yet unaddressed one-dimensional case. In addition, we pro- 
vide a critical assessment of different reciprocal-space countercharge corrections and demonstrate 
the improved accuracy of point-charge auxiliary functions in predicting the electronic levels and 
electrical response of conjugated polymers from plane-wave OD-DFT calculations. 
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I. INTRODUCTION 

Determining the electronic levels and electrical re- 
sponse of materials stands as one of the funda- 
mental limitations of density-functional theory (DFT) 
approximations.^ To address this deficiency, orbital- 
dependent density-functional theory (OD-DFT) approx- 
imations represent promising alternativess^- At present, 
the most widely used OD-DFT methods are hybrid func- 
tional approximations that consist of admixing a fraction 
of Hartree-Fock (HF) exchange into DFT functionals.3 
In hybrid functionals, orbital dependence manifests it- 
self into the nonlocal Fock contribution to the admixed 
Hamiltonian that causes the effective potential to vary 
vifith the orbital state upon which it acts. Hybrid func- 
tional methods are now well established in the field of 
electronic-structure calculations and have been shown to 
improve upon local and semilocal DFT in predicting or- 
bital properties. In parallel, self-interaction corrections 
to DFT approximations represent a second category of 
OD-DFT methods that has rapidly grown in recogni- 
tion. In the self-interaction approach, the total energy 
of the system is expressed explicitly in terms of individ- 
ual orbital densities, yielding orbital-dependent effective 
potentials. This orbital specificity is then exploited as 
an additional degree of freedom to correct unphysical 
errors inherent in the orbital-independent picture. Al- 
though progress has been gradual since the introduction 
of self-interaction corrections^ number of successful ap- 
plications have appeared recently.^ 

Despite the predictive potential of OD-DFT meth- 
ods, conventional plane-wave OD-DFT implementations 
suffer from important computational inaccuracies in 



comparison with predictions based upon atomic Slater 
functions, local Gaussian orbitals, semilocal wavelets^ 
specific discrete representations?^ and real-space finite 
differences.^ These inaccuracies, which are not intrin- 
sic to the plane-wave representation, arise from diver- 
gence singularities in the reciprocal-space summation of 
interaction contributions. Those are typically removed in 
the crudest manner by equating diverging terms to zero, 
causing finite-size errors that vanish slowly with the size 
of the computational supercelL— 

In order to cancel plane-wave errors, electrostatic cor- 
rections are necessary. For conventional DFT approxi- 
mations, many successful real-space countercharge cor- 
rections have been proposed.S^-i^ Those improve compu- 
tational accuracy at the cost of some consented increase 
in computational cost. However, in typical OD-DFT cal- 
culations, real-space countercharge corrections become 
costly due to the fact that OD-DFT functionals require 
to compute a different interaction potential, and thus a 
different correction, for each electron state at variance 
with Kohn-Sham DFT that requires instead the solution 
of one unique Hamiltonian problem per self-consistent 
iteration of the electronic-structure calculation. 

Auxiliary-function techniquesi^^— are reciprocal-space 
countercharge corrections that offer the advantage of not 
imposing any increase in computational cost in plane- 
wave OD-DFT calculations. Nevertheless, there exists 
much freedom in the possible choices of auxiliary func- 
tions and different definitions result in various levels of 
accuracy in eliminating plane- wave errors. Recently, Bro- 
qvist, Alkauskas, and Pasquarello have made an impor- 
tant step toward removing arbitrariness in the selection of 
auxiliary functions by introducing a simplified constant- 



shift correction,?^ This constant-shift correction allows 
to describe energy observables with improved accuracy. 
However, the method is not intended to rectify proper- 
ties that depend on the precise profile of the effective 
OD-DFT potential, such as the spatial confinement and 
electrical response of electron orbitals. 

In this work, we determine exact point-charge 
reciprocal-space auxiliary functions for the OD-DFT de- 
scription of molecular systems exhibiting arbitrary trans- 
lational symmetries. 

The study is organized as follows. First, we recall the 
methodological framework of reciprocal-space auxiliary- 
function corrections. Second, we examine auxiliary- 
function corrections for nonperiodic and two-dimensional 
molecular structures, thereby recovering existing correc- 
tions in simple explicit forms. We then address the un- 
treated one-dimensional case. In the last section, we 
assess the precision of the exact point-charge auxiliary- 
function correction in describing the electronic levels and 
electrical response of semiconducting oligomers of finite 
and infinite lengths. 



II. METHOD 

Plane-wave electronic-structure calculations for mate- 
rials that do not exhibit three-dimensional periodicity are 
typically carried out within the supercell approximation^ 
which consists of calculating the physical properties of an 
electronic system in a periodic finite cell and extrapolat- 
ing computed observables in the infinite-cell limit where 
interactions between the system and its artificial peri- 
odic images can be safely neglected. In this section, 
we present the framework of reciprocal-space auxiliary- 
function countercharge corrections for the cancellation of 
periodic-image errors in supercell calculations. 

In order to put matters into perspective, we consider a 
nonperiodic charge density /9od(r) confined within a su- 
percell of volume ft. Its three-dimensional periodic coun- 
terpart is written P3d(r) — J2r Pod(r-I-R) where the vec- 
tors R describe invariant translations of the superlattice. 
The Fourier components of the charge density read 

P3d(g) - ^ I drpod(r)e-'sr ^IJ drp3d(r)e-8r, (i) 

where C stands for the supercell. Within the supercell 
approximation, the potential 'y3d(r) is obtained by solv- 
ing algebraically the Poisson equation in reciprocal space, 
yielding U3d(g) = ^P3d(g) with the notable exception of 
the component ■y3d(g = 0) that is singular due to the 
diverging factor. This singularity is conventionally re- 
moved by simply omitting the singular component, i.e.. 



V3d(r) = V— P3d(g)e" 



(2) 



Physically, this crude correction can be interpreted as im- 
mersing the charge density in an artificial compensating 
jelliumj^ 



In contrast to the scheme presented above, auxiliary- 
function countercharge corrections do not rely on adding 
compensating jellium contributions to the charge density. 
They consist instead of directly transforming the periodic 
electrostatic kernel -J by inserting an auxiliary function 
4>{g) into the reciprocal-space sum: 



An 



i^3d(r) = Y.[ — + ^'^(g) P3disy 



(3) 



In Eq. ([2]), the correction A0(g) is defined in terms of 
the auxiliary function 0(g) as 



A0(g) = ^ / dr0(r)e-'S'- - ^g(g 
^ Jws 9 



(4) 



where g{g) represents the periodic Fourier decomposi- 
tion of the auxiliary density g{r) = — -^V^(/)(r) and WS 
stands for the Wigner-Seitz cell centered around the ori- 
gin. 

The auxiliary function must be chosen in order to can- 
cel divergence singularities, and in such a manner that 
V3d (r) closely reproduce the potential of the original non- 
periodic density 



wod(r) = / dr 



/Pod(r') 



(5) 



The appropriate auxiliary function can be determined by 
considering a point charge /9od(r) = '5(r). In this simple 
case, it is straightforward to show that the choice 



n 

r 
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yields the exact potential everywhere within the Wigner- 
Seitz cell. The resulting simple plane- wave correction can 
be expressed as 



A0(g) 



dr- 
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(7) 



which effectively removes the divergence singularity. 

Then, making use of the principle of superposition, 
it can be shown that the corrected potential U3d(r) is 
also exact for any given charge density pod(r) under the 
requirement that the diameter of the arbitrary density 
/3od(r) does not exceed half of that of the Wigner-Seitz 
cell, an important condition that is always tacitly as- 
sumed in the framework of reciprocal-space corrections. 

Now, the central difficulty in making use of Eq. ([7]) 
is evaluating the integral over WS due to the nonspher- 
ical WS geometry and the divergence of the integrand 
at the origin. In order to allow the evaluation of the 
reciprocal-space correction, one alternative ansatz con- 
sists of restricting the exact point-charge potential to a 
sphere of radius Tc centered at the origin)iii^i2i 



(r) 



^ if r < Tc 
if r > Tf. 



(8) 



By Fourier transform, it can be shown that Eq. 
amounts to choosing 

47r 

9 



(9) 



which also cancels out divergence singularities?^'' How- 
ever, Eq. dn]) presents the notable disadvantage of con- 
tributing oscillatory cosine terms to the potential W3d(r) 
as a result of the sharp truncation at r = Tc- Such oscil- 
latory contributions affect the convergence of calculated 
electronic observables as a function of the size of the su- 
percell and plane-wave kinetic-energy cutoff. Addition- 
ally, it should be noted that truncation schemes cannot 
be efficiently applied to elongated structures (see Sec. 
113). 

Therefore, in order to remove singularities without af- 
fecting numerical convergence, another widely used reg- 
ularization technique consists of applying a Gaussian 
smearing to Eq. ©-i^ Explicitly, the resulting Gaussian 
auxiliary function reads 



.(r) = ^erf(I 



that is, 



A0.(g) = / —erf 
/ws T 



dv 



Q 



-igr 



An 

g 



(10) 



(11) 



The Gaussian correction defined in Eq. (ITT|) is well 
behaved, correctly removes singularities in plane-wave 
sums, and does not introduce oscillatory components in 
the computation of W3d(r). Moreover, its computational 
precision can be improved systematically by decreasing 
the spread parameter a. 

Nevertheless, by examining the analytical behavior of 
Gaussian auxiliary-function corrections as a function of 
cr, Broqvist, Alkauskas, and Pasquarello have evidenced 
that the rate of convergence of the correction in the 
vanishing spread limit is poor (see Fig. 1 in Ref. 1241 ). 
Moreover, computing the limit is costly due to the very 
fine grids required in representing Gaussians of vanishing 
spread, so that, in practice, the evaluation of the exact 
correction can only be performed for a restricted num- 
ber of interpolation points. In fact, in the constant-shift 
correction of Broqvist, Alkauskas, and Pasquarello, the 
evaluation of the limit is only done at r = 0, i.e.. 



A<^(g) = 



lim 



2rj 



EAtt 



if 5 = 



if g ^ 0. 



(12) 

This simple constant-shift approach has been shown to 
improve the convergence of hybrid-functional total ener- 
gies and orbital levels as a function of supercell parame- 
ters. However, as already mentioned, the constant-shift 
auxiliary-function correction is not apt at and admittedly 
not intended to correcting electronic properties that de- 
pend on the precise profile of interaction potentials, such 



as the confinement of electronic orbitals and their re- 
sponse to an electric-field perturbation. In Sec. Illli we 
derive point-charge auxiliary-function corrections for the 
precise determination of electrostatic and exchange inter- 
action contributions. 



III. POINT-CHARGE AUXILIARY FUNCTIONS 

In this section, we present our computational proce- 
dure to obtain exact point-charge auxiliary functions for 
nonperiodic and periodic structures. The derivation con- 
sists of examining the convergence of Gaussian auxiliary- 
function corrections [Eq. (jlip ] in the vanishing spread 
limit where the correction becomes exact. ^-i 

This analysis allows to compute point-charge auxiliary 
functions at minimal computational cost. We first ap- 
ply this simplified approach to recover the non per iodic 
and two-dimensional corrections of Refs. |l^ andKLOl We 
then address the one-dimensional case, thereby comple- 
menting the derivation of point-charge auxiliary-function 
corrections for systems of arbitrary periodicity. 



A. Nonperiodic and tw^o-dimensional molecular 
structures 

We first focus on the asymptotic behavior of A0o-(r), 
which can be determined by differentiating the Fourier 
transform of Eq. dill) with respect to the squared spread, 
obtaining 






n 



,3 ^ 



(13) 



In Eq. (J13p , one can observe that the sum term vanishes 
exponentially as u goes to zero. This allows to write 
a precise modified approximation for exact point-charge 
corrections in terms of their Gaussian counterparts: 



A(/)<,=o(r) = A(/)^(r)+7ra2 



(14) 



Equation (|14p corresponds to a simple modification of 
the original Gaussian auxiliary function [Eq. (fTO|) ]: 



0.(r) 



—erf ( — 



(15) 



The simplified derivation presented above underscores 
the significance of the quadratic term ttct^ that enters 
into the expression of the original Martyna-Tuckerman 
correction. This constant contribution, which we refer 
to as the point-c harg e modification term, is only briefiy 
alluded to in Ref. 1 18 and is frequently omitted in practi- 
cal implementations of auxiliary-function schemes. ^^ To 
illustrate the importance of the point-charge modifica- 
tion term, we compare the convergence of the modified 
auxiliary function [Eq. (jlSp ] to that of the unmodified 
Gaussian auxiliary function [Eq. ((TU)) ] in Fig. [TJ demon- 
strating the considerably improved convergence of the 



modified correction to the exact point-charge limit for a 
test charge pod(r) — 6{r) in a cubic ceU of miit volume. 
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FIG. 1: Energy correction for a point charge in a cubic cell 
of volume fl = 1 a.u. with and without point-charge modi- 
fication term as a function of the spread parameter a. The 
auxiliary-function correction with point-charge modification 
converges rapidly to the exact Madelung energy of a point 
charge immersed in a jellium ^ — 2.837297 Ha. 

We now turn to the two-dimensional case. Similarly to 
Eq. p5|) . the expression of the two-dimensional correction 
reads 



<PAr) 



^'P2d,a{r), 



(16) 



where (/92d,CT(r) denotes the potential of a two-dimensional 
array of Gaussian charges that can be expressed by lon- 
gitudinal Fourier transform as 



<<52d,cr(r) 



Yl '^2d,, 



■(r±;g||)e 



'Sll"-!! 



(17) 



In Eq. ([T7)) . the vectors r|| and r^ stand for the longitu- 
dinal and transverse components of the position vector 
' " . The Fourier components ip2d.a{'<^±',S\\) sat- 



r = r| 

isfy the equation 



9r^ 



^-.gjl ) ¥'2d,a(r±;g||) = --^5id,a(r_L)e" 



(18) 

where gid.aiic±) denotes a one-dimensional Gaussian dis- 
tribution of spread a and S stands for the longitudinal 
surface area of the supercell. Now, one can determine 
the physical solution (/?2d,cr(r±; g|| = 0) of Eq. (^5)) to be 

/ \ 27r / „ /?'_L\ cr ^--i 

'P2d,a[r±;g\\ =0) = -— I r_Lcrf (^— j + -^e ^ 

(19) 
Additionally, in the general case where (7j| > 0, the 
solution of Eq. (|25l) is also obtained straightforwardly 
through a Gaussian convolution involving the point- 
charge solution 



'P2d,<T=o(r±;g||) = - 



27re" 



9\\ 



(20) 



In explicit terms, the component </32d,o-(r±; g||) can be 
written as 



(P2d.a{rr,g\\) = -j^ (2cosh(.g||rj 

Sg\\ 



e-sr^erf 



r 



^9\\ 



G 



e^r^erf^^-^ 



(21) 



As a final result, the modified Gaussian auxiliary function 
for two-dimensional molecular structures reads 



^^(r) = TTcr^ — ( r^erf ( -^ ) + -j=t 



^ Y^ ip2d,a{r±;g\\)e 






(22) 



To assess the performance of the two-dimensional 
point-charge auxiliary function, we calculate the work 
function of a model pentacene thin film consisting of 
one pentacene molecule in an orthorhombic supercell of 
fixed transverse dimensions ai = 5 A and 02 = 10 
A, and varying longitudinal dimension 03. Here, we 
use the Perdew-Zungcr (PZ) OD-DFT self-interaction 
correction"* and we employ norm-conserving local-spin 
density (LSD) pseudopotentials^ with a plane-wave en- 
ergy cutoff of 40 Ry. (Further details on the applica- 
tion of point-charge auxiliary-function corrections to pe- 
riodic systems are presented in Sec. II VI ) In Fig. [21 we 
observe that the convergence of the work function with- 
out auxiliary-function correction is very poor. Applying 
the correction is found to improve convergence; expect- 
edly, while corrected predictions are still unreliable for 
cell sizes comparable to thickness of the pentacene film, 
optimal accuracy (within a few meV) is achieved above 
30 A, that is, approximately twice the thickness of the 
film, thereby demonstrating the performance of the two- 
dimensional auxiliary-function correction. 
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FIG. 2: PZ work function of a model pentacene film as a func- 
tion of the longitudinal supercell parameter with and without 
point-charge auxiliary-function correction. 



B. One-dimensional molecular structures 

In this section, we examine molecular structures ex- 
hibiting one-dimensional translational symmetries. In 
the one-dimensional case, the modified Gaussian auxil- 
iary function can be obtained from a simple generaliza- 
tion of Eq. ([T5|) . i.e., 



(/)cr(r) = TTCr^ + ri</7id,<T(r), 



(23) 



where (/3id,cr(r) is the electrostatic potential of a one- 
dimensional periodic array of Gaussian charges. As pre- 
viously, the only potential difficulty in determining the 
point-charge auxiliary function pertains to the compu- 
tation of iy9id,CT(r). The computational procedure is ex- 
plained below. 

The one-dimensional Gaussian potential can be calcu- 
lated from the Fourier ansatz 



V'id,.(r)=^¥.id,.(ri;g||)e'«ii'-Ii 



(24) 



Indeed, substituting Eq. ([M)) into the Poisson electro- 
static problem for a one-dimensional array of Gaussian 
charges, we obtain 
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r±_ dr± dr± 



g\\ vid,^(rj_;g||) 






(25) 



which involves the 



transverse Gaussian function 

52d,cr(r±) = -^e^^^ and the distance L between 
periodic Gaussians. For the component corresponding 
to the vanishing longitudinal wavevector g|| = 0, the 
physically admissible solution readsi^ 



Ei - 



ipid,a{r±;g\\ =0) = - (-In 

(26) 
In Eq. ((26|) . Ei denotes the exponential-integral function 
and 7 stands for the Euler constant. 

We can now focus on the more difficult case where 
gii > for which it is convenient to introduce the rescaled 
radial coordinate x = g\\rj_ and to study the solution 
ka (x) of the differential equation 



dx^ 



dx 



ka{x) = -~27rx^g2,a{x), (27) 



which is related to the original solution iy9id,CT(rj_; g||) 
through 



¥'id,^(r_L;g||) =2fcg||^(5||r_L) 



(28) 



The solution of Eq. (I?7l) is slightly technical. Thus, for 
the sake of clarity, we defer the analytical derivation and 
numerical calculation of k^ (x) to Appendix \^ 



Having determined the physical solution of Eq. (|27l) . 
the final expression of the point-charge auxiliary function 
reads 



Mr) = 7ra^ + ^(-\n('l 



+ 



2n 
T 
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kg^^a{g\\rA 



)e- 
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(29) 



which completes the explicit derivation of auxiliary func- 
tions for systems of arbitrary periodicity. 

The performance of the one-dimensional point-charge 
auxiliary function is discussed and demonstrated in 
Sec.lTVl 



IV. OLIGOMERS AND POLYMERS 

In this section, we evaluate the precision of point- 
charge auxiliary-function corrections in describing the 
electronic properties of periodic and nonperiodic molec- 
ular structures. 

At the implementation stage, we have found prac- 
tical to collect point-charge auxiliary-function subrou- 
tines into a self-contained libafcc (library of auxiliary- 
function countercharge corrections) module that takes as 
an input information about the geometry, periodicity, 
and grid resolution of the supercell and returns the cor- 
rection calculated from Eqs. P3|) . (|22p . and (l^^l on the 
real-space grid. The calculated correction is then added 
to the Coulomb kernel for the computation of electro- 
static and exchange electron interactions, pseudopoten- 
tial terms, and interatomic forces. 

An early implementation of the point-charge auxiliary- 
function counterch arg e correction has been employed by 
the authors in Ref. [2^ for the OD-DFT description of iso- 
lated molecules and clusters. Here, we present a critical 
assessment of different auxiliary-function techniques in 
determining the electronic properties of extended molec- 
ular chains. Besides their known technological relevance 
to organic optoelectronics,^^ semiconducting polymers 
represent a critical test in assessing the predictive abil- 
ity of electronic-structure methodsi^"— (In fact, con- 
ventional DFT approximations systematically destabi- 
lize donor charge-carrier levels in conjugated polymers 
and strongly overestimate their longitudinal electrical re- 
sponse and optical cross section. )^i^i^ 

We consider an isolated CeHg molecule in a geometry 
recurrently studied in the literature, which consists of a 
conjugated carbon backbone with alternating simple and 
double bonds of 1.451 A and 1.339 A, respectively.2Si21i2i 

We first focus on the supercell convergence of the ion- 
ization potential that is computed as the opposite energy 
of the highest occupied molecular level. In carrying out 
these first calculations, we employ norm-conserving LSD 
pseudopotentials with a plane-wave cutoff of 50 Ry for 
the Fourier expansion of the wave functions. Electronic 
optimization is achieved using Car-Parrinello fictitious 



damped dynamics, as implemented in the CP code of the 
QUANTUM-ESPRESSO distributionJii 



size errors in comparison to the screened correction above 
9 A. 
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FIG. 3: PZ ionization potential of CsHg as a function of 
the distance between periodic replicas with and without 
reciprocal-space corrections. 



FIG. 4: PZ static longitudinal polarizability per monomer of 
CfiHg as a function of the distance between periodic replicas 
with and without reciprocal-space corrections. 



Figure [3] depicts the performance of auxiliary-function 
corrections in predicting the PZ ionization potential as a 
function of the periodic distance between artificial repli- 
cas. In Fig. [3l the truncated correction is that defined in 
Eq. (|n]), the shifted correction is identical to the constant- 
shift correction of Eq. (|12p . and the screened correc- 
tion is the auxiliary-function already available in the PW 
code of QUANTUM-ESPRESSO (version 4.3.1), which we 
have adapted to the CP code. Explicitly, this screened 
auxiliary-function correction reads 



A</>(g) 



27rCT2 if 5 = 

2 2 

A0^(g)e~^^ 



if 9 ^ 0, 



(30) 



where A0(j(g) is defined in Eq. ([TT|) and the term 2na^ 
can be obtained following the point-charge modification 
procedure described in Sec. IIIII 

The results reported in Fig. [3] confirm the slow conver- 
gence of the PZ ionization potential without correction. 
At a vacuum separation of 15 A, that is, three times 
larger than the length of CgHg, the predicted ioniza- 
tion potential is still underestimated by more than 2 eV. 
Applying the truncated correction is found to improve 
convergence substantially. However, as already men- 
tioned above, the performance of the truncated method 
for extended systems is strongly limited by the fact 
that the spherical truncation diameter should not ex- 
ceed the transverse size of the supercell, which is itself 
much smaller than the actual longitudinal extent of the 
oligomer. The constant-shift correction results in a bet- 
ter cancellation of finite-size errors, effectively reducing 
the lack of convergence to a few tenths of an eV beyond 
a distance of 9 A. In the same range, the screened cor- 
rection method is found to be remarkably precise with 
errors on the order of 0.02 meV relative to the estimated 
ionization energy. Optimal convergence is obtained with 
the point-charge correction [Eq. ((T5|) ] that halves finite- 



We now examine the static linear polarizability of 
CgHg within PZ. In these calculations, the longitudinal 
polarizability is evaluated by finite difference from the 
static dipole moment computed at a field of 0.005 a.u. 
using Berry-phase techniques i^""— The plane- wave cutoff 
energy and the oligomer geometry are kept unchanged. 
In Fig. m we observe that the truncated correction ex- 
hibits large numerical instabilities and that the constant- 
shift correction does not improve predicted polarizabili- 
ties relative to uncorrected calculations. The latter ob- 
servation is explained by the fact that the constant shift 
does not affect the profile of the orbital-dependent poten- 
tial. In contrast, both the point-charge and screened cor- 
rections improve electrical-response predictions. A closer 
comparison between the two methods turns to be in fa- 
vor of the point-charge correction in the range 3-5 A and 
in favor of the screened correction in the range 9-15 A, 
thereby illustrating the potential benefit of screening the 
point-charge auxiliary-function correction to describe the 
electrical response of molecular systems. Nevertheless, as 
discussed above (Fig. [3]) , such improvement comes at the 
cost of altering the precision of predicted orbital energies 
so that the point-charge correction should be preferred 
in general. 

To complement this comparison, we confront plane- 
wave corrected polarizability calculations to electrical 
predictions based upon local orbitals, namely, Slater-type 
orbitals (STOs) and Gaussian- type orbitals (GTOs). For 
the purpose of this analysis, it is necessary to raise the 
plane-wave energy cutoff to 60 Ry; with this compu- 
tational parameter, we verify that the longitudinal po- 
larizability of CgHg is converged within less than one 
tenth of an atomic unit. Polarizabilities are again cal- 
culated with a Berry-phase electric field of 0.005 a.u. 
using the point-charge countercharge method. The re- 
sults reported in Table U illustrate the very good per- 
formance of the point-charge correction in reproducing 



TABLE I: LSD and HP static longitudinal polarizability per 
monomer of oligoacetylenes as a function of the number of 
monomer with and without point-charge auxiliary-function 
correction, as compared with Slater-type-orbital (STO) and 
Gaussian-type-orbital (GTO) calculations. 
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Plane waves 














LSD 60 Ry 


33.5 


43.6 


57.3 


72.6 


88.8 


106.0 


HP 60 Ry 


33.1 


42.3 


52.8 


63.0 


72.3 


80.7 


Local orbitals 














LSD" STO 


32 


42 


56 


71 


87 


105 


HP*' GTO 6-31G 


37.4 


47.3 


57.2 


66.4 


74.7 


82.2 



"Reference 271. 
'Reference [34. 



Slater-type-orbital (STO) predictions. Nevertheless, it 
must be noted that there exist large deviations with 
respect to Gaussian-type-orbital (GTO) 6-31G Hartree- 
Fock predictions. A careful analysis of basis-set conver- 
gence reveals that the observed discrepancy is most likely 
due to the absence of polarization and diffuse functions 
in 6-31G. This interpretation is in line with the findings 
of Refs. |3^ and |40 and is clearly corroborated by the 
calculations presented below. 

In order to elucidate the origin of the discrepancy be- 
tween plane-wave and local-orbital predictions, we now 
examine the static polarizability of finite dimerized hy- 
drogen chains with alternating intramolecular and inter- 
molecular distances of 2 a.u. and 3 a.u., respectively. 
The electrical response of dimerized hydrogen chains is 
the subject of a wide body of literature^i"— i^i^ on as- 
sessing the relative performance of electronic-structure 
methods and describing optical saturations^ in semicon- 
ducting polymers. 
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PIG. 5: HP static longitudinal polarizability per dimer of 
dimerized hydrogen chains with and without point-charge 
auxiliary-function correction, as compared with Gaussian- 
type-orbital (GTO) 6-311G(d,p) calculations. 

Based upon GTO calculations, we can examine the 
absolute precision of our corrected plane-wave method 



in describing the electrical response of extended semi- 
conducting oligomers. Our calculations employ the same 
parameters as those used in the calculations of Table ID 
The vacuum separation between periodic chains is now of 
60 atomic units. Figure [5] demonstrates the beneficial im- 
pact of the auxiliary- function method in predicting polar- 
izabilities on a par with 6-311G(d,p) calculations within 
HF — and, by extension, within any type of hybrid func- 
tional admixture. 

To confirm this trend, we compare plane-wave LSD, 
HF, and PZ polarizabilities with GTO calculations in 
Table [Til Here, the comparison of HF results is par- 
ticularly enlightening; while 6-31G(d,p) underestimates 
the static longitudinal polarizability of dimerized hydro- 
gen chains, predictions based upon 6-311G(d,p), which 
accounts for contributions from polarization functions, 
and 6-311G-|--|-(d,p), which incorporates the effect of ad- 
ditional diffuse functions, are found to be in very close 
agreement with corrected plane-wave calculations. The 
same level of agreement is obtained for LSD and PZ pre- 
dictions with errors as low as 0.1-0.2 a.u. per dimer. 
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PIG. 6: HP static longitudinal polarizability per dimer of an 
infinite dimerized hydrogen chain (2 dimers per supercell) as 
a function of the distance between periodic replicas with non- 
periodic and one-dimensional point-charge auxiliary function 
corrections. 

We conclude with the OD-DFT description of polar- 
izability saturation in semiconducting polymers;^ The 
converged OD-DFT description of the properties of infi- 
nite chains entails to treat nonperiodic orbital-dependent 
contributions to the effective potential separately from 
periodic orbital-independent contributions. In other 
words, the problem of determining the electronic states 
of periodic structures within OD-DFT methods is akin to 
describing a localized charge density confined in a neutral 
periodic host. 

Explicitly, the electrostatic correction that we employ 
here (the one-dimensional auxiliary-function method) 
consists of computing contributions from the total 
periodic density using the one-dimensional correction 
[Eq. (j29l) ] and contributions from individual orbitals us- 
ing the nonperiodic correction [Eq. (J15p ]. This approach 
differs from the nonperiodic method that consists in- 



TABLE II: LSD, HF, and PZ static longitudinal polarizability per dimer of dimerized hydrogen chains as a function of the 
number of dimers with and without point-charge auxiliary-function correction, as compared with Gaussian-type-orbital (GTO) 
calculations at different levels of basis-set refinement. Periodic calculations are performed with 14 hydrogen units per supercell 
with a vacuum separation of ~8 A. Extrapolated GTO results are given in parentheses. 
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60 Ry 


12.4 


18.7 
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39.9 


55.4 


HF 


60 Ry 


11.9 


16.0 


18.8 


20.7 




22.0 


23.0 


23.8 


24.4 


28.7 


PZ 


60 Ry 


11.7 


16.4 


19.8 


22.2 




23.9 


25.2 


26.2 


27.0 


32.6 



LSD 
HF 



PZ 



Gaussians 

6-31H-l-G(d,p)" 

6-31G(d,p)'' 

6-311G(d,p)'' 

6-311-H-HG(d,p)" 

6-311-h-l-G(d,p)" 



11.3 
12.3 



18.8 
15.3 
16.0 
16.1 
16.5 



24.3 
18.1 
18.7 
18.9 
19.9 



28.8 
20.1 
20.6 
20.8 
22.3 



21.5 
21.9 



35.4 
22.5 
22.9 
23.1 
25.3 



37.6 
23.3 
23.6 
23.9 
26.3 



23.9 
24.2 



(28.5) 
(28.6) 



"Reference 321. 
'Reference l4ll. 



stead of computing periodic and nonperiodic contribu- 
tions indiscriminately, using the same nonperiodic auxil- 
iary function. 

Figure [H] compares the performance of the one- 
dimensional point-charge correction to that of the non- 
periodic correction. Uncorrected results based upon the 
supercell approximation of Refs. '30' and '4^ are also re- 
ported. Here, the HF electrical response of infinite hydro- 
gen chains is determined keeping calculation parameters 
unchanged, with 2 hydrogen dimers per supercell to min- 
imize computational cost. In Fig. [SJ we observe that the 
nonperiodic and one-dimensional corrections exhibit sim- 
ilar convergence behaviors. Nevertheless, we also observe 
deviations of more than 2 a.u. between the two methods 
due to the inappropriate treatment of one-dimensional 
electrostatic and exchange contributions within the non- 
periodic approach. 

The accuracy of one-dimensional calculations is further 
confirmed by direct comparison with extrapolated GTO 
calculations, as reported in the last column of Table [Hi 
As a matter of fact, plane- wave Berry-phase predictions 
are found to be in very close agreement with extrapo- 
lated 6-311G(d,p) results^ with deviations as low as 0.1 
a.u., thereby establishing the very good precision of the 
point-charge auxiliary-function correction in describing 
extended oligomers and infinite polymers without extrap- 
olation procedures. 



tions to achieve accurate convergence of electrostatic and 
exchange interactions at minimal computational cost. 

Because OD-DFT functionals require to compute a dif- 
ferent interaction potential for each electron state at vari- 
ance with DFT methods, any increase in the cost of eval- 
uating interaction terms represents a potentially critical 
computational bottleneck, limiting the interest of elabo- 
rate real-space countercharge corrections. 

To derive point-charge auxiliary functions, we have 
examined the convergence of Gaussian auxiliary func- 
tions to the exact point-charge limit. In the process, we 
have highlighted the significance of corrective contribu- 
tions that are frequently omitted in practical implemen- 
tations. We have also demonstrated the performance of 
the method in describing the frontier levels and linear po- 
larizability of molecular structures. Additionally, point- 
charge countercharge methods have been shown to put 
plane-wave OD-DFT calculations on a par with refined 
local-orbital calculations and to allow the description of, 
e.g., infinite conjugated polymers without resorting to 
delicate extrapolation procedures. 

In an effort to facilitate the incorporation of the point- 
charge auxiliary-function correction into conventional 
plane-wave codes, the correction has been implemented 
as a self-contained libafcc module.''^ In this form, we 
expect the method to prove useful in exploring further 
the accuracy of OD-DFT approximations for the descrip- 
tion of periodic molecular structures. 



V. CONCLUSION 

In summary, we have presented a reciprocal-space com- 
putational method that relies solely on plane-wave tech- 
niques and allows to study nonperiodic and periodic 
molecular structures within the predictive framework 
of OD-DFT approximations. The approach employs 
reciprocal-space point-charge auxiliary-function correc- 
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Appendix A: Derivation and computation of ka{x) 

In this appendix, we derive the analytical expression 
and explain the numerical calculation of the solution 
ka{x) of Eq. ((27)) . The function ka{x) is involved into 
the longitudinal Fourier expansion of the electrostatic po- 
tential </>id,o-(r) generated by a one-dimensional periodic 
array of Gaussian charges [Eq. ([M])]. 

Focusing first on the limit a -> where the Gaus- 
sian function transforms into a Dirac distribution, the 
solution of Eq. (I27p with vanishing asymptotic conditions 
reads 



kQ{x) = Ko{x), 



(Al) 



where Ko{x) is the modified Bessel function of the sec- 
ond kind. Now, making use of the superposition princi- 
ple, the solutions ka^(){x) can be obtained by Gaussian 
convolution: 

ko^ix) = f d('V^o(|x - y\)92dMy)- (A2) 



In practice, evaluating the above two-dimensional inte- 
gral is difficult, except in the specific case where x = 0. 
Making use of cylindrical coordinates, we can evaluate 
the function at x = to be 



kc.{0) = ~-e-Ki(- — 



(A3) 



Now, the general solution of Eq. (PT]) can be computed 
from the conventional series expansion*^ 



+ 00 



ka{x) = ^{an + bn \n{x)) a;2". 



(A4) 



n=0 



Substituting Eq. (jA4|) into Eq. ([27t yields recursive rela- 
tions of the form 



bn 
n 



fln-l + Cn-1 



bn 



bn-1 



4n^ 



{n> 1), 



(A5) 



where the coefficients c„ 



2(-l)" 

Hn-iy. 



and Co = are those 



entering into the expansion of the source term appearing 
on the right hand side of Eq. ([27]). Hence, the recursion 
defined by Eq. (jA5j) allows to express ka{x) from the 
knowledge of the first terms of the sequence qq = k^ (0) 
and 6o = whose expression is obtained from Eq. (jA3p . 



The dependence of ka{x) on x is depicted in Fig. [T] 
We observe that the values computed from the series ex- 
pansion [Eq. (jA4p ] vanish progressively upon increasing 
X until reaching a critical point where double-precision 
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FIG. 7: ka{x) computed by power series summation 
[Eq. (|A4|| ] and direct Runge-Kutta integration [Eq. (f27)l ]. 



computations become inaccurate. In this region of nu- 
merical inaccuracy, we instead resort to iterative integra- 
tion techniques for the solution of Eq. ^7} based upon 
the known asymptotic behavior 



"-a yx) — ^c 



(A6) 



of the function in the limit x -^ +00.^''' In Eq. (IA6|) . Aa 
is a constant that depends only on a and that can be 
determined straightforwardly from the requirement that 
ka{x) should not diverge at the origin. This completes 
the computation of k^ {x) . 
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still be written as reciprocal-space corrections [Eq. Q]. few tenths of an a.u., as supported by the comparison with 

It is important to note that the reciprocal-space point- refined GTO results in Sec lIVI 



